function [XS,bIS,dS] = smoothDomainBoundary( X, dx, bI, distFunct )
% SMOOTHDOMAINBOUNDARY Smooths a domain boundary for doing OMA
%
% This function takes a domain boundary and smoothes it a bit so that all
% edge segments have the same length.  This is needed for doing OMA on a
% domain - segments can't be too long and there can't be any identical
% points.
%
% The output domain boundary will have edges that are approximately the
% desired size.  In general though, the edges will not all have the same
% length because the distance along the boundary is calculated along the
% original boundary which wiggles.  Also, the length along the original
% boundary of each edge will be adjusted so that a whole number of such
% segments fits along each section of the boundary.  In rare cases, repeated
% iterations of this function might be necessary to produce the desired
% level of smoothing.
%
% Usage: [XS,bIS,dS] = smoothDomainBoundary( X, dx, bI, distfunct )
%
% Inputs:
% ------
% X = two column matrix of domain boundary coordinates.
% dx = the desired spacing among edge segments. If domain edge is in
%      Lon,Lat coordinates, dx should be in METERS!, otherwise dx should
%      be in same units as domain edge itself.
% bI = a list of indices of rows in X indicating the boundaries between open
%      and closed parts of the boundary (as generated by makeDomainBoundary
%      for example).  This list should not include the first or last point
%      of the domain boundary.  These boundaries will be respected in the
%      final smoothed domain boundary.  Defaults to [].
% distFunct = string name of function to use for calculating distance
%             between consecutive vertices.  This can be 'dist' to use the
%             Euclidean distance or 'm_idist' for Lon,Lat coordinates (m_map
%             must be on matlab path).  Defaults to 'm_idist'.
%
% Outputs:
% -------
% XS = smoothed boundary
% bIS = boundaries between open and closed sections of smoothed boundary
% dS = length of each segment of the smoothed boundary
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 	$Id: smoothDomainBoundary.m 439 2007-06-12 16:37:54Z dmk $	
%
% Copyright (C) 2007 David M. Kaplan
% License: GPL (Gnu Public License)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist( 'distFunct', 'var' )
  distFunct = 'm_idist';
end

if ~exist( 'bI', 'var' )
  bI = [];
end

% Add first and last grid point to list of edge sections
bI = [ bI(:)', size(X,1) ];

% Calculate distances
d = feval( distFunct, X(1:end-1,1), X(1:end-1,2), X(2:end,1), X(2:end,2) );
d = d(:);

[XS,bIS] = deal([]);

% Loop over pieces of open boundary
for k = 1:numel(bI)-1
  bIS = [ bIS, size(XS,1)+1 ];
  x = X(bI(k):bI(k+1),:);
  c = cumsum( [ 0; d(bI(k):bI(k+1)-1) ] );

  % Remove duplicates
  [c,I] = unique( c );
  x = x(I,:);
  
  n = max( 2, floor(c(end)/dx)+1 );
  cc = linspace( 0, c(end), n );
  
  xx = interp1( c, x, cc );
  
  XS = [ XS; xx(1:end-1,:) ];
end

XS(end+1,:) = X(end,:);

if nargout > 2
  dS = feval( distFunct, XS(1:end-1,1), XS(1:end-1,2), XS(2:end,1), ...
              XS(2:end,2) );
end
